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Progress is reported on several questions that bedevil un- 
derstanding of granular systems: (i) are the stress equations 
elliptic, parabolic or hyperbolic? (ii) how can the often- 
observed force chains be predicted from a first-principles con- 
tinuous theory? (iii) How to relate insight from isostatic sys- 
tems to general packings? Explicit equations are derived for 
the stress components in two dimensions including the de- 
pendence on the local structure. The equations are shown 
to be hyperbolic and their general solutions, as well as the 
Green function, are found. It is shown that the solutions 
give rise to force chains and the explicit dependence of the 
force chains trajectories and magnitudes on the local geome- 
try is predicted. Direct experimental tests of the predictions 
are proposed. Finally, a framework is proposed to relate the 
analysis to non-isostatic and more realistic granular assem- 
blies. 
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Granular systems have become a subject of intensive 
research in recent years both due to their enormous tech- 
nological importance and the fundamental theoretical 
challenges that they pose [1]. In particular, stress trans- 
mission has focused much attention following experimen- 
tal [2] [3] and numerical [4] observations that arching 
effects give rise to nonuniform stress fields [5] and in par- 
ticular to chain-like regions of large forces which can- 
not be straightforwardly described by conventional ap- 
proaches [6]. It has been recognized that to fully un- 
derstand this phenomenon in general granular packings 
it is essential to first understand stress transmission in 
isostatic systems [5]. Isostatic states are configurations 
of grains in which the intergranular contact forces can 
be determined directly from statics, namely, force and 
torque balance, without reference to stress-strain rela- 
tions. These states are characterized by low mean coor- 
dination number, which depends on the dimensionality 
and the roughness of the grains. These states have been 
shown to be easy to approach experimentally [7]. Sev- 
eral empirical [8] and statistical [5] [9] models have been 
proposed for the macroscopic equations that govern the 
stress field in such systems. Very recently, however, the 
two-dimensional case has been solved from first principles 
on the scale of a few grains [10]. The main result of the 
new theory is an equation that relates directly between 
the stress tensor a and a rank-two symmetric fabric ten- 
sor P which characterizes the microstructurc: 

Pxx^yy Pyy^xx ^Pxy^xy • (1) 

This 'constitutive' relation is a local manifestation of the 
torque balance condition beyond the global requirement 
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FIG. 1. The local geometry around a grain g. The vec- 
tors r*® connect contact points clockwise around each grain 
g and form loops, e.g. I, around physical voids. The vectors 
Iv^ extend from grain centers to loop centers. The geom- 
etry around the grain is characterized by the fabric tensor 
= ^jf*^-R'®, whose antisymmetric part gives the area 
= a}^ that is the shaded area. The symmetric part, P 
couples to the stress in eq. (1). 

that u = (j^ [10]. The tensor P can be defined at the 
grain level as 

V^, = \Y.{f'^1^^1^'') (2) 

where the indices z,j denote the Cartesian components 
X, the vectors r*^ and i?'^ are shown in fig. 1, and the 
sum runs over the loops / that surround grain g. Eq. (1) 
together with the conventional force and torque balance 
conditions 

V-a = 5 ; a = a^ (3) 

give a closed set of equations for the stress tensor. Here 
9 = i9x{^T 9y{^) is a position-dependent external force 
field. The only problem with eq. (1) was that the volume 
averages of Pij vanish and therefore that it couples be- 
tween the stress field and fluctuations in local geometric 
properties. This made its coarse-graining to macroscopic 
scales a non-trivial task, but this difficulty was eventually 
resolved in [11]. 

Relation (1) gives a fundamental closure of the stress 
equations and it is therefore useful to test its implications 
on two central questions: (a) Are the stress equations el- 
liptic, parabolic or hyperbohc [2] [8] [12] [13]. This ques- 
tion is fundamentally significant because elasticity theory 
predicts elliptic equations and any other form implies a 



significant departure from this paradigm, (b) How can 
the continuum theory predict and quantify the emergence 
of the experimentally observed force chains? Both these 
issues arc resolved here using the clear geometric inter- 
pretation of eq. (1) from the grain level up. This Letter 
is structured as follows: First, explicit equations for the 
stress components (7^ are derived and their form is de- 
termined. Second, the general solution of the equations 
is found and the Green function for an infinite medium is 
presented. Third, the solution is analysed and it is shown 
to give rise to force chains whose individual trajectories 
can be predicted. Next, experiments are suggested which 
can directly test the new results. Finally, a short discus- 
sion is presented on the relevance and possible extension 
of these results to general non-isostatic systems. 

To determine the stress field let us consider an isostatic 
system in mechanical equilibrium, whose microstructure 
is fully known, and solve eqs. (1) and (3). To solve for 
<7,j.y. for example, substitute ayy from eq. (1) into the 
balance eqs. (3) then eliminate axx from the two balance 
conditions by differentiating one with respect to x and 
the other with respect to y. Since the coarse-grained pij 
are a measure of the fluctuations of the geometric prop- 
erties, their gradients are neglected relative to the field 
gradients. This leaves the equations valid for disordered 
systems; in ordered periodic structures pij = indenti- 
cally. A similar manipulation for axx 

and ayy gives 

ijPxx^xx ^ ^Pxy^xy ^ Pyy^yy) ^ij — y) (^) 

where da stands for d/da and 

fxy — Vyy^yflx ~t- Pxx^xQy 

fxx — {Pxx^x ^" '^Pxydy^ 9x Pxx^yQy 

fyy iPyy^y ^ '^Pxy^^x) 9y Pyy^xQx 

are functions of the external loading. Note that these 
functions vanish identically if g is constant (e.g. constant 
gravity). It can be shown that p^y must be finite and, 
for generality, let us take Pxx and Pyy to be finite too. 
When either of these vanishes the solution to eqs. (1) 
and (3) simplifies but the qualitative behavior remains 
unchanged. A key observation is that under the following 
change of variables 




where S = a/— detP, eq. (4) takes the form 

{puu ^vv) ^ij ~ fij ' (6) 

Eqs. (6) make it possible to resolve the ongoing dispute 
regarding the nature of the stress field: They involve 
only second derivatives and so rule the parabolic form 
flat out. This suggests that diffusion-like interpretations, 
which have been proposed to explain the meandering of 
force chains, are irrelevant. This leaves two possibilities, 



either the equations are hyperbolic or they are elliptic. 
From (5) and (6) we see that the answer hinges on the 
sign of detP. When it is negative (positive) v is purely 
real (imaginary) and the equation is hyperbolic (ellip- 
tic). Now note that on the granular level P^ = P'^^ 
where P's = ^ (r^sR^s + Risfig^ jg the contribution of 

loop I to P^, and the sum runs over the loops around 

grain g (see fig. 1). It can be readily verified that (- 
detP'^) = (a'^)^ > where a'^ is the area enclosed 
by the quadrilateral shown shaded in fig. 1. However, 
detP^ ^ J2i detP'^ and therefore it does not enjoy such 
a convenient interpretation. Rather, a little algebra leads 
to 



S= HAaf -Y^{fa y,f'g).i^Rig y^Ri'9) (7) 

where = a'f . A careful consideration of the term 
under the square root leads to the conclusion that its av- 
erage over a sufficiently large area must be positive def- 
inite. For example, in a deformed honeycomb structure 
with mean separation c between centers of neighboring 
grains the average of S is 9c^/4 > 0. Nevertheless, it 
is unclear whether S can fluctuate to negative values on 
small scales. Observations made on topologically equiva- 
lent two-dimensional structures [14], such as liquid crys- 
talline foams [15] and emulsions [16] reveal no regions 
with imaginary v down to the cellular (granular) scale. 
But this may be a consequence of the isotropy of those 
systems and it is yet unclear whether small clusetrs of 
grains or foam vertices with S" < can exist in anisotropic 
systems. That said, the above result on the average of S 
means that even if such anisotropic systems exist, clus- 
ters of 'elliptic defects' cannot survive above a certain 
lengthscale and therefore eqs. (6) are indeed hyperbolic 
on macroscopic scales. 

Turning to the analysis of eqs. (6), their general solu- 
tion is 

a = i+$(?7 = v-u)+ i"*(C = v + u)+ B+ri + E'C, 

+ irfF{rj',C)dv'dC (8) 

Here -F is a rank-two matrix whose components are fij, 




and = pxy ± S. In solution (8) $(r?), *(() and the 
coefficients b^y are determined by the boundary data. 
The lines of constant 

V=\Pxxy-a.'^x]/S ; C=\pxxy-a.~x]/S (9) 

are the characteristic curves of the equations and play a 
significant role in the behavior of the stress. The Green 
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function of eq. (6) when fij is replaced by a (5-function 
at {uqjVq) within an infinite medium is 

G{u,v;uo,Vo) = ^[H{v - vq - [u - u^)) 

+ - 1^0 + (u ^ Mo))] (10) 

where H is the Heavyside step function. 

For illustration, consider a system occupying the half- 
plane a; > that is loaded along the y-axis with 
the boundary conditions ajj(x = 0,y) = Uij{y) and 
dx<^ijix — 0,y) = Vij{y). The choice of these bound- 
ary data follows the hyperbolic nature of the cqs. (6). 
The system is also presumed to be under a constant field 
g. It is convenient to first convert the boundary condi- 
tions to the u — V plane, where they become the Cauchy 
data. The solution for the stress then becomes 
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(11) 



where Ulj is the derivative of Uij with respect to its ar- 
gument. The lack of symmetry in (11) upon interchang- 
ing X and J/ is a result of the asymmetric choice of the 
variables u and v with respect to x and y. The acid 
test of this solution is that it can explain the emergence 
of force chains so often observed in planar systems [2]- 
[4]: Suppose that the boundary loading on the granular 
system is localized, namely, Vij and U-j are very nar- 
row for all aij. This is typically what happens when 
grains, which cannot form a perfect straight line at the 
boundary, are compressed by a flat surface, which local- 
izes the loading on protruding grains. From (11) we sec 
that the localized load propagates into the system along 
the characteristic curves r] = and ( = C(^, as shown 
schematically in fig. 2. In the x — y plane the charac- 
teristic curves meander due to the local dependence on 
the fluctuating geometric tensor P and the propagation 
of the forces into the system will correspondingly deviate 
from straight trajectories. For a concrete example, sup- 
pose that for the xx component the boundary conditions 



are Uxx{0,y) = e 



with d the size of the loading 



region, and 14x(0,t/) = 0. Then 
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(12) 



comprises of two bell-shaped peaks that propagate along 
the two characteristic curves. It is now possible to cal- 
culate the force along the characteristics: = a ■ tr^ 
and Fq — a ■ t(^, where f,, and are the unit tangents 
along the curves. Thresholding and visualizing now the 
force field magnitude there will emerge two lines trac- 
ing the characteristics, whose thicknesses depend on the 
threshold level. These are the force chainsl One may 
argue that force chains can be observed on scales of one 



or two grains while eqs. (1) and (3) are continuous and 
coarse-grained and so is their solution. But the only con- 
stitutive data in the solution involves the fabric tensor 
P and this tensor is well-defined down to the granular 
scale. Therefore there is no reason that solutions (11) 
should not apply also down to this scale. This analy- 
sis not only explains the force chains, it also predicts 
their individual trajectories as a function of the local 
geometry. For e xample, 1^ = (p^^,, Qi+)/v^pg^ -h (a+)^, 

k = iPxx,a~)/\/plx + ("")^- 

More general forms of boundary loading can be re- 
garded as a superposition of many localized forces such 
as the one described above. These will initiate pairs of 
chains that will propagate into the material [3]. 

When = or pyy = the characteristic curves 
take a slightly different form. For example, for p^x = 

T] = y - {1+ Pyy)x/{2pxy), C = 2/ - (1 - Pyy)x/{2pxy) 

and the functions fij change slightly. Nevertheless, the 

qualitative behavior of propagation from the boundaries 
into the system along the characteristic curves remains 
exactly the same. 

The above predictions concerning the force chains can 
be used to test the theory experimentally as follows: Ap- 
ply a localized load on the boundary of a granular (or cel- 
lular) system where it is possible to trace the force chains. 
This can be achieved, e.g., using photoelastic grains [3] 
or liquid-crystalline foams [15] [17] between crossed po- 
larizers. By scanning the local microstructure one can 
compute the fabric tensor P and then, using expressions 
(9) for 77 and C, one can compute the trajectories of the 
characteristic curves in the vicinity of the concentrated 
load. These trajectories can be checked for coincidence 
with the observed trajectories of the force chains to yield 
a straightforward test of the theory. 

Finally, what is the relevance of these results to gen- 
eral granular packings? The following offers a framework 
for the extension of the above theory to all dry granular 
systems. The ideas presented are based on insight from 
a recent experiment on inertia-free growing 2D granular 
piles of rigid and rough grains [7]. Two results of that 
experiments are of major significance: One is that the 
isostatic state (IS) can be approached arbitrarily closely. 
The other is that this state can be regarded as a critical 
point where a certain lengthscale diverges. This length- 
scale is the size of a yield front which moves ahead of 
the consolidating pile wherein grains constantly shift be- 
fore they finally consolidate. Significantly, this length 
characterizes the range of rearrangement as the system 
is perturbed locally. At the IS the mean coordination 
number is three, the material is marginally rigid, grains 
are just stable and a small displacement anywhere causes 
rearrangement far away. The range of rearrangement is 
the correlation length, ^. A critical-like behaviour has 
also been observed in other experiments [18]. As more 
contacts are made the system moves away from the IS, 
the density increases, grains are better supported and 
rearrangement is more confined. Since rearrangement is 



3 



mediated by forces it follows that the correlation length 
also describes the typical length of force chains. Thus, at 
the IS force chain are expected to span distances compa- 
rable to the system size. Indeed, infinitely long chains is 
exactly the prediction of eqs. (10) and (11). Away from 
the IS elastic domains form in which the stress equations 
arc elliptic. When a force chain is incident on an elastic 
region it can be shown that it splits into several weaker 
forces. The idea of the IS as a critical point has already 
been entertained in the literature [19] but within elas- 
ticity theory, thus failing to take into consideration the 
hyperbolic nature of force chains. It is proposed here that 
granular systems exhibiting finite force chains are neither 
at marginal rigidity nor fully elastic, but rather a mix- 
ture of these two states. This means, for example, that £, 
would decay as a power of the density difference from the 
critical density identified in [7] . Thus, for a complete the- 
ory of stresses in granular media it is essential to address 
such two-phase materials. Now that we possess a theory 
of the isostatic state, it is this author's belief that the 
extension to two-phase systems is well within reach and 
a detailed analysis along these lines is under preparation 
[20]. 

To conclude, this paper reported several results: 1) For 

isostatic systems it has been found that on large scales 
all the stress components follow an identical hyperbolic 
equation, but with different source terms. 2) The general 
solutions for the stress field were found and analysed and 
the Green function for an infinite medium was presented. 

3) It was shown that the solutions give rise to force chains 
whose individual trajectories can be predicted explicitly. 

4) Experiments were suggested which can directly test 
the validity of the new results. 5) It is proposed that 
general granular packings are in fact two-phase mixtures 
of isostatic and elastic regions and a framework to de- 
scribe such materials has been suggested. 
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